function F = demand(prc_vec,nu,a,b,par_beta)

%
% Demand for N goods given prices and parameters
%
% Strategy: determines demand by locating indifferent consumers and seeing
% whether these consumers are between [a,b] on the income distribution.
% For a given vintage, utility is increasing in income.
% 

N = length(prc_vec);
F = zeros(N,1);
rec_pts = zeros(N,3);

% what does the richest guy buy?
    [util_ptr, high_bnd] = max(nu.*(b-prc_vec));
% what does the poorest guy buy?
    [util_ptr, low_bnd] = max(nu.*(a-prc_vec));
    
% Case where everyone buys one product
if(low_bnd==high_bnd)
    F(low_bnd) = 1.0;

% 2 or more products are sold
else
    j=low_bnd;
    while j<=high_bnd-1 
            
        % for product j, calculate indifferent consumers for each
        % higher utility
        indiff = zeros(high_bnd-j,1);
        for i=j+1:high_bnd
            indiff(i-j) = (nu(i)*prc_vec(i)-nu(j)*prc_vec(j))/(nu(i)-nu(j));
        end    
        % Selecting the vintage with the lowest indifference consumer
        [inc_ptr, n_ptr] = min(indiff);
        % check: indifferent consumer's income is too low or too high
        if (inc_ptr<a || inc_ptr>b)
            disp('Problem in routine'); disp([nu prc_vec]); disp(indiff); disp([high_bnd low_bnd]);
            disp([inc_ptr,a,b,low_bnd,high_bnd,j,i]);
            pause
        else
            rec_pts(j,:) = [inc_ptr, j, j+n_ptr];   
            j = j+n_ptr;            
        end
    end
    % converting rec_pts into demand
    % need to transform rec_pts on [0,1] range, for beta function
    indx = (rec_pts(:,1)>0);
    rec_ptsT = indx.*(rec_pts(:,1)-a)/(b-a);
    
    indx = find(rec_pts(:,1)>0);
    ll_indx = length(indx);
    
    % only two vintages are sold
    if(ll_indx==1)
        tmp = cdf('beta',rec_ptsT(indx,1),par_beta(1),par_beta(2));
        F(rec_pts(indx,2)) = tmp*1.0;
        F(rec_pts(indx,3)) = (1-tmp)*1.0;
    else
        for j=1:length(indx)
            if(j==1)
                F(rec_pts(indx(j),2)) = cdf('beta',rec_ptsT(indx(j),1),par_beta(1),par_beta(2))*1.0;                 
            elseif(j==length(indx))
                tmp0 = cdf('beta',rec_ptsT(indx(j),1),par_beta(1),par_beta(2));
                tmp1 = cdf('beta',rec_ptsT(indx(j-1),1),par_beta(1),par_beta(2));
                F(rec_pts(indx(j),2)) = (tmp0-tmp1)*1.0;
                F(rec_pts(indx(j),3)) = (1-tmp0)*1.0;
            else
                tmp0 = cdf('beta',rec_ptsT(indx(j),1),par_beta(1),par_beta(2));
                tmp1 = cdf('beta',rec_ptsT(indx(j-1),1),par_beta(1),par_beta(2));
                F(rec_pts(indx(j),2)) = (tmp0-tmp1)*1.0;
            end
        end
    end

end    
       
    % check on demand
    if (abs(sum(F)-1.0)>0.00001)
        disp('Demand does not sum to mrkt_size (1.0)');
        disp([sum(F) 1.0]); F, prc_vec,pause
    end
    





            
